using Plots, Distributions, DynamicGridsConway’s Game of Life
Game of life (code in order of presentation)
What are we building?
We are implementing GOL, resulting in something like the following
@time out = sim!(ArrayOutput(rand(Bernoulli(0.3), (100,100)), tspan=1:100), Life());
anim = @animate for i in 1:100
heatmap(out[i], c=:binary, aspectratio=1)
end;
gif(anim, fps=10) 2.568484 seconds (6.82 M allocations: 374.577 MiB, 4.54% gc time, 99.55% compilation time)
┌ Info: Saved animation to
│ fn = /home/michael/phd/GOL/tmp.gif
└ @ Plots /home/michael/.julia/packages/Plots/1KWPG/src/animation.jl:114
Starting by defining the rules of GOL
Rules
- Any living cell with < 2 neighbors dies.
- Any living cell with two or three live neighours lives
- Any cell with > 3 neighbors dies
- Any dead cell with 3 neighbors becomes alive
We want to design it not as a one-off simulation, by like we are designing a software package for GOL.
Why do this? It leads to better, more reliable, and more reusable code
The Board struct
struct Board{T<:Number}
grid::Matrix{T}
end
Base.size(b::Board) = size(b.grid)
Base.getindex(b::Board, ci::CartesianIndex) = b.grid[ci]
Base.setindex!(b::Board, val, ci::CartesianIndex) = setindex!(b.grid, val, ci)
Base.eachindex(board::Board) = CartesianIndices(size(board.grid))
Base.similar(b::Board{T}) where T = Board(zeros(T, size(b)))Base.in(i::CartesianIndex, board::Board) = begin
x,y = size(board)
i[1] <= x && i[1] > 0 && i[2] > 0 && i[2] <= y
endBoard(p = 0.3, sz = (100,100)) = Board(rand(Bernoulli(p), sz))Board
The simulation loop
function simulate(init::Board, numtimesteps=100)
timeseries = [similar(init) for _ in 1:numtimesteps]
timeseries[begin] = init
for t in 2:numtimesteps
update!(timeseries[t-1], timeseries[t])
end
return timeseries
endsimulate (generic function with 2 methods)
So we’re done? No, we need to implement the update! function
function update!(oldboard, newboard)
for cellindex in eachindex(oldboard)
update!(oldboard, newboard, cellindex)
end
endupdate! (generic function with 1 method)
and now the update! function for each index
function update!(oldboard, newboard, cellindex)
st = oldboard[cellindex]
numneighbors = countneighbors(oldboard, cellindex)
newboard[cellindex] = nextstate(st, numneighbors)
endupdate! (generic function with 2 methods)
Now we need to final functions: countneighbors and nextstate.
function countneighbors(board, index)
offsets = filter(x->x!=CartesianIndex(0,0), CartesianIndices((-1:1,-1:1)))
neighborindices = [index + o for o in offsets]
filter!(x-> x ∈ board, neighborindices)
return sum(board.grid[neighborindices])
endcountneighbors (generic function with 1 method)
function nextstate(state, numlivingneighbors)
if !state
return numlivingneighbors == 3
elseif numlivingneighbors == 2 || numlivingneighbors == 3
return true
elseif numlivingneighbors < 2 || numlivingneighbors > 3
return false
end
endnextstate (generic function with 1 method)
@time simulate(Board(0.1,(50,50)), 100); 0.128976 seconds (990.10 k allocations: 121.109 MiB, 26.04% gc time)
@time simulate(Board(0.1,(50,50)), 100); 0.126039 seconds (990.10 k allocations: 121.109 MiB, 28.53% gc time)
function animate(timeseries)
pltsettings = (
aspectratio=1,
cbar=:none,
frame=:box,
c=:binary)
return @animate for t in 1:length(timeseries)
heatmap(timeseries[t].grid;
pltsettings...,
xlim=(1,size(timeseries[begin])[1]),
ylim=(1,size(timeseries[begin])[2])
)
end
end animate (generic function with 1 method)
seq = simulate(Board(0.2,(50,50)), 100);
gif(animate(seq), fps=10)┌ Info: Saved animation to
│ fn = /home/michael/phd/GOL/tmp.gif
└ @ Plots /home/michael/.julia/packages/Plots/1KWPG/src/animation.jl:114
Adding Patterns
abstract type Pattern end
Base.size(s::Pattern) = s.sz
grid(s::Pattern) = s.shapegrid (generic function with 1 method)
struct Glider <: Pattern
sz
grid
end
function Glider()
mat = [false false true;
true false true;
false true true];
return Glider(size(mat), mat)
end
(g::Glider)() = g.gridstruct Blinker <: Pattern
sz
grid
end
function Blinker()
mat = [false true false;
false true false;
false true false]
return Glider(size(mat), mat)
end
(b::Blinker)() = b.gridstruct Pulsar <: Pattern
sz
grid
end
(p::Pulsar)() = p.grid
function Pulsar()
mat = [0 0 0 0 1 0 0 0 0 0 1 0 0 0 0;
0 0 0 0 1 0 0 0 0 0 1 0 0 0 0;
0 0 0 0 1 1 0 0 0 1 1 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 1 1 0 0 1 1 0 1 1 0 0 1 1 1;
0 0 1 0 1 0 1 0 1 0 1 0 1 0 0;
0 0 0 0 1 1 0 0 0 1 1 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 1 0 0 0 1 1 0 0 0 0;
0 0 1 0 1 0 1 0 1 0 1 0 1 0 0;
1 1 1 0 0 1 1 0 1 1 0 0 1 1 1;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 1 0 0 0 1 1 0 0 0 0;
0 0 0 0 1 0 0 0 0 0 1 0 0 0 0;
0 0 0 0 1 0 0 0 0 0 1 0 0 0 0]
return Pulsar(size(mat), mat)
end Pulsar
# Note: no boundscheckfunction add!(board::Board, shape::T, base::Tuple) where T<:Pattern
sz = size(shape)
upper = base .+ sz .- 1
board.grid[base[1]:upper[1],base[2]:upper[2]] .= shape()
board
endadd! (generic function with 1 method)
add!(Board(), Glider(), (10,10))Board{Bool}(Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 1 … 1 0; 0 1 … 0 0])
Base.:+(board::Board, s::Shape) = begin
add!(board, s)
return board
end board = Board(zeros(Bool, 70,70))
add!(board, Glider(), (10,10))
add!(board, Blinker(), (35,50))
add!(board, Pulsar(), (30,30))Board{Bool}(Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0])
heatmap(board.grid, c=:binary, aspectratio=1);out = simulate(board, 500);gif(animate(out), fps=10)┌ Info: Saved animation to
│ fn = /home/michael/phd/GOL/tmp.gif
└ @ Plots /home/michael/.julia/packages/Plots/1KWPG/src/animation.jl:114
DynamicGrids.jl
GOL
init = rand(Bernoulli(0.3), (100,100));@time out = sim!(ArrayOutput(init, tspan=1:100), Life()); 0.024749 seconds (561 allocations: 1.036 MiB, 72.43% compilation time)
const DEAD, ALIVE, BURNING = 1, 2, 3
neighbors_rule = let prob_combustion=0.0001, prob_regrowth=0.01
Neighbors(Moore(1)) do data, neighborhood, cell, I
if cell == ALIVE
if BURNING in neighborhood
BURNING
else
rand() <= prob_combustion ? BURNING : ALIVE
end
elseif cell == BURNING
DEAD
else
rand() <= prob_regrowth ? ALIVE : DEAD
end
end
end[31mNeighbors[39m(
f = var"#9#10"{Float64, Float64}
neighborhood = Moore{1, 2, 8, Nothing}
)
init = fill(ALIVE, 200, 200)
output = ArrayOutput(init;
tspan=1:500,
fps=15,
);@time sim!(output, neighbors_rule); 1.356233 seconds (3.26 M allocations: 177.713 MiB, 11.26% gc time, 73.44% compilation time)
cs = cgrad([:black, :green, :orange]);anim = @animate for i in 1:length(output)
heatmap(output[i], c=cs, clim=(1,3))
end Animation("/tmp/jl_NczJPE", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png" … "000491.png", "000492.png", "000493.png", "000494.png", "000495.png", "000496.png", "000497.png", "000498.png", "000499.png", "000500.png"])
gif(anim)┌ Info: Saved animation to
│ fn = /home/michael/phd/GOL/tmp.gif
└ @ Plots /home/michael/.julia/packages/Plots/1KWPG/src/animation.jl:114